Análisis del paisaje

Autores/as

Afiliaciones

Víctor Gauto

GISTAQ (UTN-FRRe)

IIDTHH (UNNE, CONICET)

Instituto Gulich (UNC, CONAE)

Sofía Paná

Instituto Gulich (UNC, CONAE)

El paquete usado es {rflsgen}, Neutral Landscape Generator with Targets on Landscape Indices.

Precaución

Antes que nada hay que ejecutar los siguientes comandos:

options(
  java.parameters = c("-XX:+UseConcMarkSweepGC", "-Xmx36g")
)
gc()

Esto permite aumentar la memoria disponible para las funciones de {rflsgen}.

Los paquetes usados son los siguientes:

library(tidyverse)
library(gt)
library(landscapemetrics)
library(NLMR)
library(rflsgen)
library(terra)

1 Lectura de datos

Inicialmente hay que leer los ráster de uso de suelo (r) y modelo digital de elevación (dem).

r <- rast("datos/r.tif")
dem <- rast("datos/dem.tif")

Las clases de r son:

Valor	Tipo
0	No data
10	Forest
20	Scrubland
30	Grassland
40	Agriculture
50	Urban
60	Sparse vegetation - bare soil
80	Water body

Agrego las categorías a r e incorporo colores. Visualizo ambos ráster.

d <- read_delim(
  file = "datos/LULCcode.txt",
  delim = "$ ",
  col_names = c("value"),
  show_col_types = FALSE
) |>
  separate_wider_delim(
    delim = " ",
    cols = value,
    names = c("value", "tipo"),
    too_many = "merge"
  ) |>
  mutate(value = as.integer(value)) |>
  filter(value != 0) |>
  mutate(col = c(
    "#E41A1C", "#377EB8", "#4DAF4A", "#984EA3", "#FF7F00", "#A65628", "#F781BF"
  ))

col_df <- d |>
  select(-tipo) |>
  as.data.frame()

# incorporo categorías y colores
levels(r) <- d
coltab(r) <- col_df

2 Extracción de la estructura del paisaje

Se seleccionan las clases de interés del ráster de cobertura del suelo.

clases <- c(30, 40, 50) # Grassland, Agriculture, Urban

estructura <- flsgen_extract_structure_from_raster(
  raster_file = r,
  focal_classes = clases,
  connectivity = 8
)

A partir de estructura y dem se genera un paisaje, indicando dependencia del terreno, conectividad, rugosidad y parámetros de posición.

Para visualizar el efecto de la conectividad se generan 4 escenarios con valores crecientes: 0.1, 0.3, 0.7 y 1.

f_conectividad <- function(x) {
  print(paste0("Conectividad: ", x))
  
  r_land <- flsgen_generate(
    estructura,
    terrain_file = dem,
    terrain_dependency = x,
    connectivity = 8,
    epsg = "EPSG:22174",
    resolution_x = res(r)[1],
    resolution_y = res(r)[2],
    x = ext(dem)$"xmin",
    y = ext(dem)$"ymax"
  )
  
  writeRaster(r_land, filename = paste0("rasters/land-", x, ".tif"))
}

conectividad <- c(.1, .3, .7, 1)
walk(conectividad, f_conectividad)

Los escenarios generados son:

El paquete {rflsgen} genera múltiples escenarios pero manteniendo constantes las métricas del paisaje.

Para el cálculo de las métricas se utiliza el paquete {landscapemetrics}. Se calcularon: Mean of patch area y Number of patches.

f_area_mn <- function(x) {
  lsm_c_area_mn(land_raster[[x]], directions = 8) |>
    mutate(conectividad = conectividad[x]) |> 
    mutate(metrica = "Mean of patch area")
}

df_area_mn <- map(1:length(conectividad), f_area_mn) |> 
  list_rbind()

write_tsv(df_area_mn, "indices/df_area_mn.tsv")

f_np <- function(x) {
  lsm_c_np(land_raster[[x]], directions = 8) |>
    mutate(conectividad = conectividad[x]) |> 
    mutate(metrica = "Number of patches")
}

df_np <- map(1:length(conectividad), f_np) |> 
  list_rbind()

write_tsv(df_np, "indices/df_np.tsv")

Los valores se obtienen para cada clase elegida previamente. Se observa que para cada valor de conectividad las métricas coinciden.

read_tsv("indices/df_area_mn.tsv", show_col_types = FALSE) |> 
  filter(class != -1) |> 
  mutate(value = round(value, 2)) |> 
  inner_join(d_r, by = join_by(class == value)) |> 
  pivot_wider(
    names_from = conectividad,
    values_from = value,
    names_prefix = "Conect: "
  ) |> 
  select(Tipo = tipo, starts_with("conect")) |> 
  gt() |> 
  tab_header("Mean of patch area") |> 
  tab_style(
    locations = cells_body(columns = starts_with("Conect")),
    style = cell_text(font = "JetBrains Mono", color = "black")
  ) |> 
  fmt_number(
    columns = everything(),
    dec_mark = ",",
    decimals = 2,
    sep_mark = "."
  )

Mean of patch area
Tipo	Conect: 0.1	Conect: 0.3	Conect: 0.7	Conect: 1
Grassland	3,52	3,52	3,52	3,52
Agriculture	29,54	29,54	29,54	29,54
Urban	1,88	1,88	1,88	1,88

read_tsv("indices/df_np.tsv", show_col_types = FALSE) |> 
  filter(class != -1) |> 
  mutate(value = round(value, 2)) |> 
  inner_join(d_r, by = join_by(class == value)) |> 
  pivot_wider(
    names_from = conectividad,
    values_from = value,
    names_prefix = "Conect: "
  ) |> 
  select(Tipo = tipo, starts_with("conect")) |> 
  gt() |> 
  tab_header("Number of patches") |> 
  tab_style(
    locations = cells_body(columns = starts_with("Conect")),
    style = cell_text(font = "JetBrains Mono", color = "black")
  ) |> 
  fmt_number(
    columns = everything(),
    dec_mark = ",",
    decimals = 0,
    sep_mark = "."
  )

Number of patches
Tipo	Conect: 0.1	Conect: 0.3	Conect: 0.7	Conect: 1
Grassland	45.871	45.871	45.871	45.871
Agriculture	12.308	12.308	12.308	12.308
Urban	15.311	15.311	15.311	15.311

En la generación de los escenarios pueden suceder dos errores:

• Falta de memoria, al no contar con la suficiente capacidad de procesamiento. Para solucionarlo se puede reducir el tamaño de la región de interés, aumentar el tamaño de píxel o disminuir el número de clases.

• Incapacidad de generar el ráster, mostrando el siguiente mensaje:

Could not generate a raster satisfying the input landscape structure

Esto puede solucionarse eligiendo únicamente las clases de interés.

3 Estructura del paisaje propia

Genero el caso I y el caso II, cada uno con dos clases y métricas del paisaje diferentes.

Luego, para verificar los escenarios generados, calculo las métricas y verifico que se cumplan las condiciones dadas.

Defino dos clases y sus métricas.

tibble(
  Clase = c("Clase 1", "Clase 2"),
  `Number of patches` = c("[700, 800]", "[800, 1100]"),
  `Patch area` = c("[5500, 14600]", "[9500, 11060]"),
  `Proportion of landscape` = c("10%", "-")
) |> 
  gt() |> 
  cols_label(
    `Number of patches` = md("Number of<br>patches"),
    `Proportion of landscape` = md("Proportion of<br>landscape"),
    "Clase" = ""
  ) |> 
  tab_style(
    locations = cells_body(columns = -Clase),
    style = cell_text(font = "JetBrains Mono", color = "black")
  ) |> 
  tab_header(title = md("**Caso I**"))

Caso I
	Number of patches	Patch area	Proportion of landscape
Clase 1	[700, 800]	[5500, 14600]	10%
Clase 2	[800, 1100]	[9500, 11060]	-

tibble(
  Clase = c("Clase A", "Clase B"),
  `Number of patches` = c("[700, 800]", "[700, 800]"),
  `Patch area` = c("[15500, 114600]", "[15500, 114600]"),
  `Proportion of landscape` = c("-", "25%")
) |> 
  gt() |> 
  cols_label(
    `Number of patches` = md("Number of<br>patches"),
    `Proportion of landscape` = md("Proportion of<br>landscape"),
    "Clase" = ""
  ) |> 
  tab_style(
    locations = cells_body(columns = -Clase),
    style = cell_text(font = "JetBrains Mono", color = "black")
  ) |> 
  tab_header(title = md("**Caso II**"))

Caso II
	Number of patches	Patch area	Proportion of landscape
Clase A	[700, 800]	[15500, 114600]	-
Clase B	[700, 800]	[15500, 114600]	25%

# CASO I
cls_1 <- flsgen_create_class_targets(
  class_name = "Clase 1",
  NP = c(700, 800),
  AREA = c(5500, 14600),
  PLAND = c(10, 10)
)

cls_2 <- flsgen_create_class_targets(
  class_name = "Clase 2",
  NP = c(800, 1100),
  AREA = c(9500, 11060)
)

# CASO II
cls_A <- flsgen_create_class_targets(
  class_name = "Clase A",
  NP = c(700, 800),
  AREA = c(15500, 114600)
)

cls_B <- flsgen_create_class_targets(
  class_name = "Clase B",
  NP = c(700, 800),
  AREA = c(15500, 114600)
)

Las métricas y los valores tienen que ser coherentes y no pueden adoptar cualquier número. En caso de indicar valores incorrectos, se muestra el siguiente mensaje de error:

User targets cannot be satisfied

Creo un conjunto con las clases y obtengo la estructura, utilizando el dem como máscara.

# CASO I
objetivos <- flsgen_create_landscape_targets(
  mask_raster = dem,
  classes = list(cls_1, cls_2)
)

estructura_custom <- flsgen_structure(objetivos)

# CASO II
objetivos2 <- flsgen_create_landscape_targets(
  mask_raster = dem,
  classes = list(cls_A, cls_B)
)

estructura_custom2 <- flsgen_structure(objetivos2)

Genero un ráster a partir de las clases creadas y su estructura ubicados sobre el dem original.

r_custom <- flsgen_generate(
  estructura_custom,
  terrain_file = dem,
  terrain_dependency = 1,
  connectivity = 8,
  epsg = "EPSG:22174",
  resolution_x = res(r)[1],
  resolution_y = res(r)[2],
  x = ext(dem)$"xmin",
  y = ext(dem)$"ymax"
)

writeRaster(r_custom, "rasters/r_custom.tif", overwrite = TRUE)

r_custom2 <- flsgen_generate(
  estructura_custom2,
  terrain_file = dem,
  terrain_dependency = 1,
  connectivity = 8,
  epsg = "EPSG:22174",
  resolution_x = res(r)[1],
  resolution_y = res(r)[2],
  x = ext(dem)$"xmin",
  y = ext(dem)$"ymax"
)

writeRaster(r_custom2, "rasters/r_custom2.tif", overwrite = TRUE)

Visualizo los resultados de ambos casos.

r_custom <- rast("rasters/r_custom.tif")

d_custom <- tibble(
  value = as.integer(c(0, 1)),
  tipo = c("Clase 1", "Clase 2"),
  col = c("#E41A1C", "#4DAF4A")
) |> 
  add_row(
    value = as.integer(-1),
    tipo = "X",
    col = "grey95"
  )

col_df_custom <- d_custom |>
  select(-tipo) |>
  as.data.frame()

# incorporo categorías y colores
levels(r_custom) <- d_custom
coltab(r_custom) <- col_df_custom

plot(r_custom, axes = FALSE, mar = c(1, 1, 1, 4), main = "Caso I")
r_custom2 <- rast("rasters/r_custom2.tif")

d_custom2 <- tibble(
  value = as.integer(c(0, 1)),
  tipo = c("Clase A", "Clase B"),
  col = c("#E41A1C", "#4DAF4A")
) |> 
  add_row(
    value = as.integer(-1),
    tipo = "X",
    col = "grey95"
  )

col_df_custom2 <- d_custom2 |>
  select(-tipo) |>
  as.data.frame()

# incorporo categorías y colores
levels(r_custom2) <- d_custom2
coltab(r_custom2) <- col_df_custom2

plot(r_custom2, axes = FALSE, mar = c(1, 1, 1, 4), main = "Caso II")

Verifico que el ráster creado tenga las métricas previamente definidas para cada clase.

read_tsv("datos/df_custom.tsv", show_col_types = FALSE) |> 
  filter(class != -1) |> 
  mutate(
    clase = if_else(class == 0, "Clase 1", "Clase 2")
  ) |> 
  pivot_wider(
    names_from = metrica,
    values_from = value,
    id_cols = clase
  ) |> 
  select(clase, `Number of patches`, `Percentage of landscape of class`) |> 
  gt() |> 
  cols_label(
    `Number of patches` = md("Number of<br>patches"),
    `Percentage of landscape of class` = md("Percentage of<br>landscape of class"),
    "clase" = ""
  ) |> 
  tab_style(
    locations = cells_body(columns = -clase),
    style = cell_text(font = "JetBrains Mono", color = "black")
  ) |> 
  fmt_number(
    columns = 2,
    dec_mark = ",",
    decimals = 0,
    sep_mark = "."
  ) |> 
  fmt_number(
    columns = 3,
    dec_mark = ",",
    decimals = 1,
    sep_mark = "."
  ) |> 
  tab_header(title = md("**Caso I**"))

Caso I
	Number of patches	Percentage of landscape of class
Clase 1	791	10,0
Clase 2	806	11,5

read_tsv("datos/df_custom2.tsv", show_col_types = FALSE) |> 
  filter(class != -1) |> 
  mutate(
    clase = if_else(class == 0, "Clase A", "Clase B")
  ) |> 
  pivot_wider(
    names_from = metrica,
    values_from = value,
    id_cols = clase
  ) |> 
  select(clase, `Number of patches`, `Percentage of landscape of class`) |> 
  gt() |> 
  cols_label(
    `Number of patches` = md("Number of<br>patches"),
    `Percentage of landscape of class` = md("Percentage of<br>landscape of class"),
    "clase" = ""
  ) |> 
  tab_style(
    locations = cells_body(columns = -clase),
    style = cell_text(font = "JetBrains Mono", color = "black")
  ) |> 
  fmt_number(
    columns = 2,
    dec_mark = ",",
    decimals = 0,
    sep_mark = "."
  ) |> 
  fmt_number(
    columns = 3,
    dec_mark = ",",
    decimals = 1,
    sep_mark = "."
  ) |> 
  tab_header(title = md("**Caso II**"))

Caso II
	Number of patches	Percentage of landscape of class
Clase A	701	16,3
Clase B	753	25,0